# data vizlibrary(ggplot2)library(ggh4x)library(ggOceanMaps)library(patchwork)library(viridis)library(ggpubr)library(grid)library(viridis)library(ggnewscale)library(gridExtra)library(png)# tablelibrary(gt)library(gtable)# maplibrary(ggmap)library(ggsn)library(sf)library(sp)library(smoothr)# kernel densitylibrary(eks)library(ggdensity)# statlibrary(Hmisc)library(circular)library(CircStats)# charlibrary(stringr)# solar angle and timezonelibrary(suncalc)library(lutz)# data wranglinglibrary(tidyr)library(dplyr)library(data.table)library(purrr)library(forcats)library(lubridate)
2 Setting up custom function
2.1windrose
Code
windrose <-function(data_to_plot,grid =NULL,set_title =NULL,legend_position ="none",facet_wrap = F,facet_grid = F,hour_sunset =20,hour_sunrise =7) {# this code comes from Roxanne uniqhours <-1:24* (360/24)# trick to align hours om the graph data_to_plot <-rbind(data_to_plot[2:nrow(data_to_plot), ], data_to_plot[1, ])for (i in1:24) {# turn hours to radiansif (i ==1) { temp <-rep(uniqhours[i], data_to_plot$nb_ind_hour[i]) day_night <-rep("night", data_to_plot$nb_ind_hour[i]) } else { temp <-c(temp, rep(uniqhours[i], data_to_plot$nb_ind_hour[i])) day_night <-c(day_night, rep(if_else(between(i, hour_sunrise, hour_sunset), "day", "night"), data_to_plot$nb_ind_hour[i] )) } } data2 <-data.frame(direction = temp) deg <-15# choose bin size (degrees/bin) dir.breaks <-seq(0- (deg /2), 360+ (deg /2), deg) # define the range of each bin# assign each direction to a bin range dir.binned <-cut(data2$direction,breaks = dir.breaks,ordered_result =TRUE )# generate pretty lables dir.labels <-as.character(c(seq(0, 360- deg, by = deg), 0))# replace ranges with pretty bin labelslevels(dir.binned) <- dir.labels# Assign bin names to the original data set data2$dir.binned <- dir.binned# add origin if anyif (facet_wrap =="origin") { data2$origin <-unique(data_to_plot$origin) } elseif (facet_wrap =="trip") { data2$trip <-unique(data_to_plot$trip) }# add trip and origin if anyif (facet_grid) { data2$origin <-unique(data_to_plot$origin) data2$trip <-unique(data_to_plot$trip) }# set up max value maxvalue <-20# initialise the plot plt.dirrose_2 <-ggplot()# check if gridif (!is.null(grid)) { plt.dirrose_2 <- plt.dirrose_2 +geom_hline(yintercept = grid,colour ="grey20",linewidth = .2 ) } plt.dirrose_2 <- plt.dirrose_2 +geom_vline(xintercept =c(seq(1, 24, 2)),colour ="grey30",linewidth =0.2 ) +# 24 vertical lines at center of the 30? ranges.geom_hline(yintercept = maxvalue,colour ="black",linewidth = .5 ) +# Darker horizontal line as the top border (max).# On top of everything we place the histogram bars.geom_bar(data = data2,aes(x = dir.binned, fill = day_night),width =1,colour ="black",linewidth =0.3 ) +# geom_bar(data = data2, aes(x = dir.binned2), width = 1, colour="black", size = 0.3,fill="salmon",alpha=0.9) +scale_x_discrete(drop =FALSE,labels =c(0, "", 2, "", 4, "", 6, "", 8, "", 10, "", 12, "", 14, "", 16, "", 18, "", 20, "", 22, "" ) ) +scale_fill_manual(values =c("white", "darkgrey"), name ="Time of day") +labs(x ="Time (hours)", y ="Count", title = set_title) +coord_polar(start =-(deg /2) * (pi /180))# if facetif (facet_wrap =="origin") { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(. ~ origin) } elseif (facet_wrap =="trip") { plt.dirrose_2 <- plt.dirrose_2 +facet_wrap2(trip ~ .,strip.position ="right" ) }# if facetif (facet_grid) { plt.dirrose_2 <- plt.dirrose_2 +facet_grid2(trip ~ origin) }# Wraps the histogram into a windrose plt.dirrose_2 <- plt.dirrose_2 +theme_bw() +theme(legend.position = legend_position,panel.grid.minor =element_blank(),panel.grid.major =element_blank(),panel.background =element_blank(),legend.key.size =unit(0.5, "cm") )# returnreturn(plt.dirrose_2) }
2.2matlab_to_posix
Code
# function to convert matlab date into R posixmatlab_to_posix <-function(x, timez ="UTC") { days <- x -719529# 719529 = days from 1-1-0000 to 1-1-1970 secs <- days *86400# 86400 seconds in a dayreturn(as.POSIXct(strftime(as.POSIXct(secs,origin ="1970-1-1",tz ="UTC" ),format ="%Y-%m-%d %H:%M:%S",tz ="UTC",usetz =FALSE ),tz = timez ))}
3 Import data
Let’s load data_dive, i.e. the output of data_wrangling.qmd, and filter only on animals leaving from Ano Nuevo.
Code
# import the datadata_dive <-readRDS("../export/data_dive.rds")# filter on seals departing from ANNUdata_dive <- data_dive %>%filter(DepartureLocation =="ANNU")
We need to define a custom function to calculate the summarize version of each column differently, i.e. weighted mean and sum.
Code
compute_table_1_a <-function(x) {# if the sum of the column is below 100if (sum(x) <100) {# it's definitely a proportion# so we weighted average these y <-weighted.mean(x, c(1068746, 2207361)) y <-paste0(round(y *100, 1), "%")return(y)# otherwise we just sum the column } else { y <-prettyNum(sum(x), big.mark =",")return(y) }}
# custom function to summarize each column differentlycompute_table_1_b <-function(x) {# if the column is numericif (all(is.numeric(x))) {# and the sum is below 100if (sum(x) <100) {# it is time# so we "circular" weighted average them y <- psych::circadian.mean(c(rep(x[1], 237),rep(x[2], 187) ),na.rm = T ) y <-round(y, 1) } else {# otherwise it is just number# so we "just" weighted average them y <-weighted.mean(x, c(237, 187)) y <-round(y, 1) } } else { y <-"\u00b1" }return(y)}
Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.
# seals
All dives
Benthic dives
% Benthic dives per area
Total
w/ loc
w/o loc
Total
Proportion
Shark
Orca
Total
Overlap
Post Breeding
237
1,068,746
1,042,451
26,295
36,443
3.4%
79.8%
100.0%
100.0%
79.8%
Post Molting
187
2,207,361
2,032,528
174,833
44,594
2.0%
90.5%
100.0%
100.0%
90.5%
OVERALL
424
3,276,107
3,074,979
201,128
81,037
2.5%
87%
100%
100%
87%
Code
# display second tabletable_1_b
Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals
# dives
# benthic dives
Hour departure
Hour arrival
Mean
±
SD
Mean
±
SD
Mean
±
SD
Mean
±
SD
Post Breeding
4509.5
±
904.4
155.1
±
275.7
18.7
±
1.8
14.8
±
2.5
Post Molting
11804.1
±
1379.4
238.5
±
471.3
18.9
±
1.7
12.9
±
2.0
OVERALL
7726.7
±
1113.8
191.9
±
362
18.8
±
1.8
13.9
±
2.3
Notes
Table 1: Summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals.
Table 2: Individual-level summary of dataset characteristics for Post Breeding and Post Molting foraging trips in northern elephant seals
4.2 Figure 1
Individual are not chosen completely randomly, we picked randomly individuals that performed ~180 dives per 36h.
Code
# to choose an individual ~180 dives for 36 hoursdata_dive %>%# rename date as datetimemutate(datetime = date) %>%# then by individualgroup_by(id) %>%# calculatemutate( # the difference time with the first divediff_start =abs(as.numeric(difftime(first(datetime), datetime,units ="days" ))),# the difference time with the last divediff_end =abs(as.numeric(difftime(last(datetime), datetime,units ="days" ))) ) %>%filter(diff_start <= nb_days | diff_end <= nb_days) %>%mutate(phase =factor(if_else(diff_start <= nb_days, "Departure", "Arrival"),levels =c("Departure", "Arrival") )) %>%group_by(id, trip) %>%summarise(nb_dives =length(DiveNumber)) %>%arrange(nb_dives) %>%View()
Code
# import the csv exporttdr_data <-lapply(list.files(path ="../data/",pattern ="*_iknos_raw_data.csv",full.names = T ),function(x) {# extract id id_name <-last(strsplit(strsplit(x, split ="_")[[1]][1],split ="/" )[[1]])# load file dt <- readr::read_csv(x)# add id dt %>%mutate(id =as.numeric(id_name)) }) %>%bind_rows()# choose the lagnb_days <-1.5# compute the datadata_plot <- tdr_data %>%# convert time into date, a Posixmutate(datetime =matlab_to_posix(time),date =as.Date(datetime) ) %>%# add beginning and end date of the tripmerge( ., data_dive %>%group_by(id, trip) %>%summarise(date_low =min(date),date_high =max(date),.groups ="drop" ),by ="id",all.x =TRUE ) %>%# identify what to keepmutate(to_keep =if_else(between(datetime, date_low, date_high), T, F)) %>%# keep only between the beginning and the endfilter(to_keep) %>%# sort to make sure the difftime are being done correctlyarrange(id, time) %>%# merge to get lat and lonmerge( ., data_dive %>%group_by(id, date =as.Date(date)) %>%summarise(lat =median(Lat, na.rm = T),lon =median(Long, na.rm = T),.groups ="drop" ),by =c("id", "date"),all.x =TRUE ) %>%# then by individualgroup_by(id) %>%# calculatemutate( # the difference time with the first divediff_start =abs(as.numeric(difftime(first(datetime), datetime,units ="days" ))),# the difference time with the last divediff_end =abs(as.numeric(difftime(last(datetime), datetime,units ="days" ))) ) %>%filter(diff_start <= nb_days | diff_end <= nb_days) %>%mutate(phase =factor(if_else(diff_start <= nb_days, "Departure", "Arrival"),levels =c("Departure", "Arrival") ))# add bathysetDT(data_plot)setDT(data_dive)data_dive[, datetime := date]setkeyv(data_plot, c("id", "datetime"))setkeyv(data_dive, c("id", "datetime"))data_plot <-copy(data_dive) %>% .[, .(id, datetime, bathy, DiveTypeName)] %>% .[data_plot, roll = T]# get sunrise sunset informationsun_data <- data_plot %>%select(date, lat, lon, id, phase) %>%unique() %>%ungroup() %>%mutate(getSunlightTimes(data = ., keep =c("sunrise", "sunset"))) %>%# fuck it, getSunlightTimes gives the sunset of next day... so we subtract 1mutate(sunset =if_else( # if sunset is the next dayabs(as.numeric(difftime(as.POSIXct(date), sunset, units ="days") )) >1,# subtract 1 sunset - (1*60*60*24),# otherwise leave it that way sunset ))# trim if night beging before the trip or end after they arrivesun_data <-merge( sun_data %>%select(date, lat, lon, id, phase, sunrise, sunset), data_plot %>%group_by(id, phase, trip) %>%summarise(date_low =min(datetime),date_high =max(datetime),.groups ="drop" ),by =c("id", "phase"),all.x =TRUE, ) %>%mutate(# make sure sunset is not outside the 36h periodcorrect_sunset =if_else(# if notbetween(sunset, date_low, date_high), sunset,# if so, replace by one extremumif_else(sunset < date_low, date_low, date_high) ),# make sure sunrise is not outside the 36h periodcorrect_sunrise =if_else(# if notbetween(sunrise, date_low, date_high), sunrise,# if so, replace by one extremumif_else(sunrise < date_low, date_low, date_high) ) ) %>%# shift everything to start 2000-01-01 00:00:00mutate(correct_sunset_shift = correct_sunset - ( date_low -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") ),correct_sunrise_shift = correct_sunrise - ( date_low -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") ) )# plotdeparture_dp <- data_plot %>%filter(phase =="Departure") %>%group_by(id, phase) %>%# make the date starts on 2000-01-01 00:00:00mutate(datetime = datetime - (first(datetime) -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") )) %>%ungroup() %>%mutate(`Dive Type`=if_else(DiveTypeName =="Benthic", "Benthic", "Other")) %>%ggplot(aes(x = datetime, y = depth, col =`Dive Type`, group = id)) +scale_color_manual(values =c("Benthic"="#008c49", "Other"="grey10")) +geom_rect(data = sun_data %>%filter(phase =="Departure"),aes(xmin = correct_sunset_shift,xmax = correct_sunrise_shift,ymin =-Inf,ymax =+Inf ),alpha =0.4,inherit.aes =FALSE ) +# geom_path(aes(x = datetime, y = abs(bathy))) +geom_path() +scale_x_datetime(# date_breaks = "6 hours",# date_labels = "+ %Hh",breaks =seq.POSIXt(as.POSIXct("2000-01-01 06:00:00", tz ="UTC"),as.POSIXct("2000-01-02 06:00:00", tz ="UTC"),by ="6 hour" ),labels =c("+6h", "+12h", "+18h", "+24h", "+30h"),position ="top" ) +scale_y_reverse() +labs(x ="Time", y ="Depth (m)") +# facet_nested(trip + id ~ phase,facet_nested(id ~ phase,scales ="free",# independent = "x" ) +theme(strip.placement ="outside",strip.text.y =element_blank(),legend.position ="top",panel.grid.minor =element_blank(),panel.border =element_blank(),panel.background =element_blank(),panel.grid.major =element_line(color ="grey90"),axis.line.x =element_line(color ="black",# arrow = arrow(length = unit(0.2, "lines"), type = "closed") ),axis.line.y =element_line(color ="black",arrow =arrow(length =unit(0.2, "lines"),type ="closed",ends ="first" ) ),axis.title.x =element_blank() )arrival_dp <- data_plot %>%filter(phase =="Arrival") %>%group_by(id, phase) %>%# make the date starts on 2000-01-01 00:00:00mutate(datetime = datetime - (first(datetime) -as.POSIXct("2000-01-01 00:00:00", tz ="UTC") )) %>%ungroup() %>%mutate(`Dive Type`=if_else(DiveTypeName =="Benthic", "Benthic", "Other")) %>%ggplot(aes(x = datetime, y = depth, col =`Dive Type`, group = id)) +scale_color_manual(values =c("Benthic"="#008c49", "Other"="grey10")) +geom_rect(data = sun_data %>%filter(phase =="Arrival"),aes(xmin = correct_sunset_shift,xmax = correct_sunrise_shift,ymin =-Inf,ymax =+Inf ),alpha =0.4,inherit.aes =FALSE ) +# geom_path(aes(x = datetime, y = abs(bathy))) +geom_path() +scale_x_datetime(# date_breaks = "6 hours",# date_labels = "+ %Hh",breaks =seq.POSIXt(as.POSIXct("2000-01-01 06:00:00", tz ="UTC"),as.POSIXct("2000-01-02 06:00:00", tz ="UTC"),by ="6 hour" ),labels =c("-30h", "-24h", "-18h", "-12h", "-6h"),position ="top" ) +scale_y_reverse() +labs(x ="Time", y ="Depth (m)") +# facet_nested(trip + id ~ phase,facet_nested(id ~ phase,scales ="free",# independent = "x" ) +theme(strip.placement ="outside",legend.position ="top",panel.grid.minor =element_blank(),panel.border =element_blank(),panel.background =element_blank(),panel.grid.major =element_line(color ="grey90"),axis.line.x =element_line(color ="black",arrow =arrow(length =unit(0.2, "lines"), type ="closed") ),axis.line.y =element_blank(),axis.ticks.y =element_blank(),axis.text.y =element_blank(),axis.title.y =element_blank(),axis.title.x =element_blank() )fig_label_x <-ggplot(data.frame(l ="Time", x =1, y =1)) +geom_text(aes(x, y, label = l)) +theme_void() +coord_cartesian(clip ="off")
Figure 1: Snapshot of departure and arrival dive profiles for six seals. Departure and arrival were constrained to 36 hours. Shading indicates night and coloration indicates dive type.
Notes
TODO: Add three dots between Departure and Arrival
4.3 Figure 2
Let’s create a table specific for this figure that only contains for each individual:
the first dive
the last dive
all benthic dives
Code
# let's add a column with the local_timedata_windrose <- data_dive %>%# get rid of data without location informationfilter(!is.na(Lat)) %>%# then by individualgroup_by(id) %>%# keep only the first, last date, but also all benthic divesfilter(date ==min(date) | date ==max(date) | DiveTypeName =="Benthic") %>%# ungroupungroup() %>%# change colnames to make getSunlightTimes worksmutate(date =as.Date(date_tz),lat = Lat,lon = Long ) %>%# calculate hour of sunset and hour of sunrisemutate(getSunlightTimes(data = .,keep =c("sunrise", "sunset"),tz =tz_lookup_coords(lat = Lat,lon = Long,method ="accurate" ) ) ) %>%# extract hours of sunset an sunrisemutate(h_sunrise =hour(sunrise),h_sunset =hour(sunset) )
Code
hour_time_2_radian <-function(x, high =24, low =0) {# https://pingouin-stats.org/build/html/generated/pingouin.convert_angles.html#pingouin.convert_angles ptp <- high - low rad <- x * (2* pi) / ptp rad <- (rad + pi) %% (2* pi) - pireturn(rad)}
Figure 2: Circular histogram plots display the local timing (in hours) for the first and last recorded dives (respectively departure and arrival) of adult female northern elephant seals (n = 403) during their Post Breeding and Post Molting foraging trips.
Note
n = 403, because we removed individuals without location data (filter(!is.na(Lat)))
Figure 3: Spatial overlap of benthic dives performed by northern elephant seals (n = 401), and the shark (yellow polygon) and orca (red polygon) distribution area along the west coast (adapted from Jorgenson, et al. 2019) with a kernel density estimation representing the occurrence of all the benthic dives.
n = 401, because we removed individuals without location data and those that do not have any benthic dives (filter(!is.na(Lat) & DiveTypeName == "Benthic")) => it means that two individuals did not have any dives classified as benthic dives… (id = 2016014 and id = 2018002)
4.5 Figure 4
Code
# https://mycolor.space/?hex=%23008C49&sub=1# build palette colorcols_ind_1 <-c("#008c49", "#00c9d3", "#bffcfd", "#458081")# test for color blind# => https://davidmathlogic.com/colorblind/#%23008C49-%2300C9D3-%23BFFCFD-%23458081cols_ind_2 <-c("#008c49", "#002e00", "#00c9d3", "#deb887")# test for color blind# => https://davidmathlogic.com/colorblind/#%23008C49-%23002E00-%2300C9D3-%23DEB887
Code
# main plotbathy_prop <- data_dive %>%# keep only rows with location datafilter(!is.na(Lat)) %>%# keep only negative bathyfilter(bathy <0) %>%# and remove outliersfilter(bathy >-6000) %>%# positive bathymutate(bathy =abs(bathy)) %>%# create class of bathymetrymutate(bathy_class =cut( bathy,seq(0, 6000, 400),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(bathy_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per bathy_classgroup_by(bathy_class) %>%# to get the Proportion of different dive types per bathy_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = bathy_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Bathymetry (m)",y ="Dive type proportion (%)" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="top",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # second plot# bathy_hist <- data_dive %>%# # keep only negative bathy# filter(bathy < 0) %>%# # and remove outliers# filter(bathy > -6000) %>%# # create class of bathymetry# mutate(bathy_class = fct_rev(cut(# bathy,# seq(-6000, 0, 400),# ordered_result = T,# dig.lab = 4# ))) %>%# # calculate by bath_class and animal# group_by(bathy_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = bathy_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## # the actual plot# bath_plot <- bathy_hist / bathy_prop + plot_layout(heights = c(1, 10))
Code
# main plotdist_coast_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_coast >0) %>%# remove outliersfilter(dist_coast *1000*111<=1900) %>%# create class of bathymetrymutate(dist_class =cut(# convert decimal degree/1000 dist_coast *1000*111,seq(0, 1900, 100),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the Proportion of different dive types per dist_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from the coast (km)",y ="Dive type proportion" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="none",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # the second plot# dist_coast_hist <- data_dive %>%# # filter out animal without dist_coast# filter(dist_coast > 0) %>%# # remove outliers# filter(dist_coast * 1000 * 111 <= 1900) %>%# # create class of bathymetry# mutate(dist_class = cut(# # convert decimal degree/1000# dist_coast * 1000 * 111,# seq(0, 2000, 100),# ordered_result = T,# dig.lab = 4# )) %>%# # calculate by bath_class and animal# group_by(dist_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = dist_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## the actual plot# dist_coast_plot <- dist_coast_hist / dist_coast_prop + plot_layout(heights = c(1, 10))
Code
# the main plotdist_dep_prop <- data_dive %>%# filter out animal without dist_coastfilter(dist_dep >0) %>%# remove outliersfilter(dist_dep /1000<4200) %>%# create class of bathymetrymutate(dist_class =cut( dist_dep /1000,seq(0, 4200, 300),ordered_result = T,dig.lab =4 )) %>%# calculate by bath_class and animalgroup_by(dist_class, DiveTypeName) %>%# the number of divessummarise(N =n(), .groups ="drop") %>%# and divide by the total number of dives per dist_classgroup_by(dist_class) %>%# to get the Proportion of different dive types per dist_classmutate(Proportion = N /sum(N)) %>%# ungroup => not required but that let the dataset cleanungroup() %>%# then plotggplot(aes(x = dist_class, y = Proportion)) +# the areageom_area(aes(fill = DiveTypeName, group = DiveTypeName)) +# orientation of x labelsguides(x =guide_axis(angle =45)) +# format y axisscale_y_continuous(labels =function(x) {paste0(x *100, "%") } ) +labs(x ="Distance from departure (km)",y ="Dive type proportion" ) +# scale_fill_viridis(# "Dive Type",# option = "plasma",# discrete = T,# direction = -1# ) +scale_fill_manual("Dive Type",values = cols_ind_1 ) +# position legend at the toptheme(legend.position ="none",axis.title.x =element_text(margin =margin(t =15)),panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )),plot.tag.position =c(-0.03, 1) )# # the second plot# dist_dep_hist <- data_dive %>%# # filter out animal without dist_coast# filter(dist_coast > 0) %>%# # remove outliers# filter(dist_dep / 1000 < 4200) %>%# # create class of bathymetry# mutate(dist_class = cut(# dist_dep / 1000,# seq(0, 4200, 300),# ordered_result = T,# dig.lab = 4# )) %>%# # calculate by bath_class and animal# group_by(dist_class, DiveTypeName) %>%# # the number of dives# summarise(N = n(), .groups = "drop") %>%# # ggplot# ggplot(aes(# x = dist_class,# y = N,# fill = DiveTypeName,# width = 0.5# )) +# geom_bar(# stat = "identity",# position = "dodge"# ) +# scale_fill_viridis(# option = "plasma",# discrete = T,# direction = -1# ) +# theme_void() +# theme(legend.position = "top")## # the actual plot# dist_dep_plot <- dist_dep_hist / dist_dep_prop + plot_layout(heights = c(1, 10))
# display(guide_area() + bathy_prop + dist_coast_prop + dist_dep_prop) +theme(plot.margin =unit(c(0, .02, -1, .05), units ="npc")) +plot_annotation(tag_levels ="A") +plot_layout(nrow =5, guides ="collect")grid::grid.draw(grid::textGrob(ylab, x =0.07, y =0.55, rot =90))
Figure 4: Proportion of each dive type performed by northern elephant seals in relation to (A) bathymetry, (B) distance from the coast, and (C) distance from departure (defined as the first recorded location).
Notes
TODO: change the color code to match the color we attributed to benthic dives (greens?)
n = 403, because we removed individuals without location data (filter(!is.na(Lat)))
4.6 Figure 5
Code
# the last plotprop_at_sea <- data_dive %>%group_by(id) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# order by seals/datearrange(id, nb_days_departure) %>%# perform our calculation per sealgroup_by(id) %>%# calculate of the Proportion of time since departuremutate(perc_time_at_sea =round( nb_days_departure /max(nb_days_departure) *100, 1 )) %>%# filter on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of Proportion of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of Proportion of time at seagroup_by(day_class) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class ) %>%# plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="#008c49") +guides(x =guide_axis(angle =45)) +scale_y_continuous(limits =c(0, 1),# breaks=seq(0,0.6,0.1),labels =function(x) { scales::percent(abs(x), 1) } ) +annotate("text", x =5.3, y =0.4, label ="Departure") +annotate("text", x =15, y =0.5, label ="Arrival") +# geom_curve(# data = tibble(# x1 = c(4, 17),# x2 = c(1.5, 19.5),# y1 = c(0.4, 0.5),# y2 = c(0.31, 0.55)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = 0.3# ) +geom_curve(data =tibble(x1 =c(4),x2 =c(1),y1 =c(0.4),y2 =c(0.33) ),aes(x = x1,y = y1,xend = x2,yend = y2 ),arrow =arrow(length =unit(0.08, "inch")),linewidth =0.5,color ="gray20",curvature =0.3 ) +geom_curve(data =tibble(x1 =c(16),x2 =c(19.7),y1 =c(0.5),y2 =c(0.57) ),aes(x = x1,y = y1,xend = x2,yend = y2 ),arrow =arrow(length =unit(0.08, "inch")),linewidth =0.5,color ="gray20",curvature =-0.56 ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Proportion of benthic dives" )
Code
# displayprop_at_sea
Figure 5: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea. Time spent at sea is classified into intervals of five percent for graphical representation.
Notes
?
n = 422, because we removed individuals that do not have any benthic dives (i.e. we keep only those having benthic dives, filter(DiveTypeName == "Benthic")) => it means that two individuals did not have any dives classified as benthic dives… (id = 2016014 and id = 2018002)
# data use to compute kernel density estimationdf_kernel_dens <- data_dive %>%# only with location datafilter(!is.na(Lat)) %>%# select only the required columnsselect(lat = Lat, long = Long, id, DiveTypeName) %>%# create col to nicely display dive typemutate(origin =paste(DiveTypeName, "dives"))
Code
# transform data into sf objectdf_kernel_dens_sf <-st_as_sf( df_kernel_dens,coords =c("long", "lat"),crs =st_crs(4326))# make it's group_by origindf_kernel_dens_sf <-group_by(df_kernel_dens_sf, origin)# kernel density estimationdf_kernel_dens_sf_kde <-st_kde(df_kernel_dens_sf,H =diag(c( MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 1] ), MASS::bandwidth.nrd( sf::st_coordinates(df_kernel_dens_sf)[, 2] ) ) /4)^2)# https://github.com/r-spatial/sf/issues/1762sf::sf_use_s2(FALSE)# plotto_plot <- trip_zoom_out +# rename axislabs(x ="Longitude", y ="Latitude") +# new scalenew_scale_fill() +# kernelgeom_sf(data =st_get_contour(# geospatial kernel df_kernel_dens_sf_kde,# probabilitiescont =c(50, 80, 95, 99) ),# displayaes(fill =label_percent(contlabel)),alpha =0.7 ) +# same colour barscale_fill_viridis_d(option ="plasma") +# legendlabs(fill ="Probs") +# no display alphaguides(alpha ="none") +# facet by originfacet_wrap(. ~ origin)# reorder layersto_plot$layers <-if (with_bathy) { to_plot$layers[c(1, 2, 4, 3)]} else { to_plot$layers[c(1, 3, 2)]}# displayto_plot
Figure 6: Map of kernel density estimation representing the spatial distribution of Benthic, Drift, Forage, and Transit dives performed by northern elephant seals throughout their trip to sea. Probs represents different levels of probability density.
Notes
lmk!
4.7.3 Figure 2
Code
# the last plotprop_at_sea_trip <- data_dive %>%group_by(id, trip) %>%# time difference between the first date and the others (telling R to take diff b/w first date and all other dates)mutate(nb_days_departure =trunc(as.numeric(difftime( date, first(date),units ="days" )))) %>%# group by day/dive_type/sealgroup_by(nb_days_departure, DiveTypeName, id, trip) %>%# count number of divessummarise(nb_daily_dives_type =n(), .groups ="drop") %>%# group by day/sealgroup_by(nb_days_departure, id) %>%# count the total number of dives per day/sealmutate(nb_daily_dives =sum(nb_daily_dives_type)) %>%# order by seals/datearrange(id, nb_days_departure, trip) %>%# perform our calculation per sealgroup_by(id, trip) %>%# calculate of the Proportion of time since departuremutate(perc_time_at_sea =round( nb_days_departure /max(nb_days_departure) *100, 1 )) %>%# filter on benthic divesfilter(DiveTypeName =="Benthic") %>%# add class of Proportion of time at seamutate(day_class =cut(perc_time_at_sea, seq(0, 100, 5),include.lowest = T )) %>%# by class of Proportion of time at seagroup_by(day_class, trip) %>%# calculate the proportion of daily benthic divessummarise(nb_benthic_class =sum(nb_daily_dives_type),nb_total_class =sum(nb_daily_dives),prop_class = nb_benthic_class / nb_total_class ) %>%# plotggplot( .,aes(x = day_class, y = prop_class) ) +geom_bar(stat ="identity", fill ="#008c49") +guides(x =guide_axis(angle =45)) +scale_y_continuous(limits =c(0, 1),# breaks=seq(0,0.6,0.1),labels =function(x) { scales::percent(abs(x), 1) } ) +facet_wrap(. ~ trip) +# annotate("text", x = 5.3, y = 0.4, label = "Departure") +# annotate("text", x = 15, y = 0.5, label = "Arrival") +# geom_curve(# data = tibble(# x1 = c(4),# x2 = c(1),# y1 = c(0.4),# y2 = c(0.33)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = 0.3# ) +# geom_curve(# data = tibble(# x1 = c(16),# x2 = c(19.7),# y1 = c(0.5),# y2 = c(0.57)# ),# aes(# x = x1,# y = y1,# xend = x2,# yend = y2# ),# arrow = arrow(length = unit(0.08, "inch")),# linewidth = 0.5,# color = "gray20",# curvature = -0.56# ) +theme(panel.grid.major =element_blank(),panel.grid.minor =element_blank(),panel.background =element_blank(),axis.line =element_line(arrow =arrow(angle =20,length =unit(.10, "inches"),type ="closed" )) ) +labs(x ="Time spent at sea (%)",y ="Proportion of benthic dives" )
Code
# displayprop_at_sea_trip
Figure 7: Proportion of benthic dives performed during trip at sea expressed as a percentage of time spent at sea during post breeding and post molting trip. Time spent at sea is classified into intervals of five percent for graphical representation.